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Abstract 

In the present paper we consider a 2-D shallow-water equations (SWE) model on 
a /3-plane solved using an alternating direction fully implicit (ADI) finite-difference 
scheme on a rectangular domain. The scheme was shown to be unconditionally- 
stable for the linearized equations. 

The discretization yields a number of nonlinear systems of algebraic equations. 
We then use a proper orthogonal decomposition (POD) to reduce the dimension of 
the SWE model. Due to the model nonlinearities, the computational complexity 
of the reduced model still depends on the number of variables of the full shallow - 
water equations model. By employing the discrete empirical interpolation method 
(DEIM) we reduce the computational complexity of the reduced order model due 
to its depending on the nonlinear full dimension model and regain the full model 
reduction expected from the POD model. 

To emphasize the CPU gain in performance due to use of POD/DEIM, we also 
propose testing an explicit Euler finite difference scheme (EE) as an alternative to 
the ADI implicit scheme for solving the swallow water equations model. 



We then proceed to assess the efficiency of POD/DEIM as a function of num- 
ber of spatial discretization points, time steps, and POD basis functions. As was 
expected, our numerical experiments showed that the CPU time performances of 
POD/DEIM schemes are proportional to the number of mesh points. Once the num- 
ber of spatial discretization points exceeded 10000 and for 90 DEIM interpolation 
points, the CPU time was decreased by a factor of 10 in case of POD/DEIM implicit 
SWE scheme and by a factor of 15 for the POD/DEIM explicit SWE scheme in com- 
parison with the corresponding POD SWE schemes. Moreover, our numerical tests 
revealed that if the number of points selected by DEIM algorithm reached 50, the 
approximation errors due to POD/DEIM and POD reduced systems have the same 
orders of magnitude supporting the theoretical results existing in the literature. 

Keyword: shallow water equations; proper orthogonal decomposition; reduced-order 
models (ROMs); finite difference methods; discrete empirical interpolation method (DEIM); 

1 Introduction 

The shallow water equations are the simplest form of the equations of motion that can 
be used to describe the horizontal (motion) structure of the atmosphere. They describe 
the evolution of an incompressible and inviscid fluid in response to gravitational and 
rotational accelerations and their solutions represent East West propagating Rossby waves 
and inertia - gravity waves. 

To avoid the limitations imposed by the Courant Friedrichs - Lewy (CFL) stability 
conditions restricting the time steps in explicit finite difference approximations, implicit 
scheme must be considered. We propose here the alternating direction implicit (ADI) 
method introduced by Gustafsson in [T| . Linear and nonlinear versions of ADI scheme may 
be found in studies proposed by Fairweather and Navon in [2] and Navon and De Villiers in 
[3] . Kreiss and Widlund in jl] established the convergence of alternating direction implicit 
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methods for elliptic problems. Such methods reduce multidimensional problem to systems 
of one dimensional problems (Douglas and Gunn [5], Yanenko [6j and Marchuck [7]). 

The nonlinear algebraic systems corresponding to the discrete model were solved using 
the quasi-Newton method proposed in Gustafsson jilj. This quasi- Newton method per- 
forms an LU decomposition done every M-th time step, where M is a fixed integer. Since 
back substitution is a fast operation the scheme will be efficient as long as the number of 
iterations is small. 

The major issue in large scale complex modelling is that of reducing the computational 
cost while preserving numerical accuracy. Among the model reduction techniques, the 
proper orthogonal decomposition (POD) method provides an efficient means of deriving 
the reduced basis for high-dimensional nonlinear flow systems. The POD method has been 
widely and successfully applied to signal analysis and pattern recognition as Karhunen- 
Loeve, statistics as principal component analysis (PCA), geophysical fluid dynamics and 
meteorology as empirical orthogonal functions (EOF) etc. The POD method was applied 
also to SWE model and we mention here the work of Cao et al. [8], Vermeulen and 
Heemink [9], Daescu and Navon pU] and Altaf et al. [TT] . 

In this paper we reduced the dimension of the SWE model by employing the POD 
method. However due to the nonlinearities of the implicit SWE model the computational 
complexity of the reduced shallow water equations model still depends on the number 
of variables of the full shallow - water equations model. To mitigate this problem, we 
apply the discrete empirical interpolation method (DEIM ) to address the reduction of the 
nonlinear components and thus reduce the computational complexity by implementing 
the POD/DEIM method. 

DEIM is a discrete variant of the empirical interpolation method (EIM) proposed by 
Barrault et al. [12] for constructing an approximation of a non-affine parameterized func- 
tion, which was proposed in the context of reduced-basis discretization of nonlinear partial 
differential equations. The application was suggested and analysed by Chaturantabut and 
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Sorensen in [15] , [TS] , [To"] and [16J. 

The paper is organized as follows. In Section 2 we introduce the Gustafsson ADI fully 
implicit method applied to the shallow water equations model and briefly describe its 
algorithmic components since they are already available in archived literature. In Section 
3 we describe in some detail the snapshot POD procedure and its implementation to the 
ADI method for the SWE model. Section 4 addresses the snapshot POD combined with 
DEIM methodology and provides the detailed algorithmic description of the DEIM imple- 
mentation. In Section 5 we present the numerical experiments related to the POD/DEIM 
procedure for both explicit and implicit schemes applied to the SWE models. 

The POD/DEIM procedure amounts to replace orthogonal projection with an interpo- 
lation projection of the nonlinear terms that requires the evaluation of only a few selected 
components of the nonlinear terms. 

We evaluate the efficiency of DEIM as a function of number of spatial discretiza- 
tion points, time steps and basis functions for this quadratically nonlinear problem and 
additional studies about the conservation of the integral invariants of the SWE, root 
mean square errors(RMSEs) and correlation coefficients between full model, POD and 
POD/DEIM systems were performed. 

2 Brief description of the Gustafsson ADI method. 

In meteorological and oceanographic problems, one is often not interested in small time 
steps because the discretization error in time is small compared to the discretization error 
in space. The fully implicit scheme considered in this paper is first order in both time 
and space and it is stable for large CFL condition numbers (we tested the stability of the 
scheme for a CFL condition number equal to 7.188). It was proved by Gustafsson in [1] 
that the method is unconditionally stable for the linearized version of the SWE model. 
Here we shortly describe the Gustafsson shallow water alternating direction implicit 
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method (Gustafsson pQ, Fairweather and Navon [2J, Navon and de Villiers [3J). We are 
solving the SWE model using the /3-plane approximation on a rectangular domain. 



_ = J 4 W _ + BW _ + G(!/K 



< x < L, < ?/ < L>, t G (0,t/], 



(1) 



where w = (u,v,<f)) T is a vector function, w,t> are the velocity components in the x and 
y directions, respectively, h is the depth of the fluid, g is the acceleration due to gravity 
and = 2yfgh. 

The matrices A, B and C are expressed 
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f = f + P(y-D/2), = !/e[0,4 



with / and constants. 

We assume periodic solutions in the x-direction 



w(x,y,t) = w(x + L,y,t), x = 0, ye [0,D], t e (0,t f ], 
while in the y— direction 

v(x, 0, t) = v(x, D,t) — 0, a; e [0, L], t G (0, £/]. 
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Initially w(x,y,0) = ip(x,y), ip : R x R ->• R, (x,y) G [0,L] x [0,D]. Note that no 
boundary conditions are necessary for u and at y = 0, D. 

Now we introduce a mesh of iV x • N y equidistant points on [0, L] x [0, D], with Ax = 
L/(N X — 1), Ay = D/(N y — 1). We also discretize the time interval [0,tf] using NT 
equally distributed points and At = tf/(NT — 1). Next we define vectors of unknown 
variables of dimension n xy = N x ■ N y containing approximate solutions such as 

u(t n ) fau(xi,yj,t n ),v(t n ) « v(x h y jt t n ), <f>(t n ) m (j)(x i: y 5 , t n ) eR Hxy , 

i = l,2,..,N x , j = 1,2, .., N y , n = 1,2,. .NT. 

Then Gustafsson's nonlinear ADI finite difference shallow water equations scheme 
(ADI FD SWE) is defined by 

I. First step - get solution at t n+ i 




with " *" denoting the componentwise multiplication, f is a A^-dimensional vector storing 

the Coriolis components f(yj), j = 1, 2, ..N y and the nonlinear functions F n , F 12 , F 21 , F 22 , F 31 , F 32 : 

R n *y x R nxy — > R nxy are defined as follows 

F n (u, 0) = u * A x u A x <f), 
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F 12 (u, v) = v * A y u, F 2 i(u, v) = u * A x v, 



F 22 {v, <fi) = v * A y v + ^<fi* A y 4>, 
F 31 (u, cj>) = ^0 * A x u + u * A x (f>, 
F 3 2(v, (j)) = ^0 * A y v + v* A y 4), 



where A x ,A y e M" 1 ^" 19 are constant coefficient matrices for discrete first-order and 
second-order differential operators which take into account the boundary conditions. 
II. Second step - get solution at t n+ i 

u(Wi) + y Fi2 ( u (Vi),v(Wi)J - Y^izl^ilf* v (^+i) = u (*n+|)- 

N x 

^Fn^u(t n+ i),0(t n+ i)J, 
v(Wi) + Y F 2 2 (v(t n ),0(t n )) = v(t n+ i) - ^F 21 (u(t n+ p,v(t n+ i)j- 



(3) 



2 



0(*n+l) + Y F 32(v(t n+ l),0(t n+ l)^) = 0(t n+ i) " ^F 31 (u(t n+ i),0(t n+ i 

The nonlinear systems of algebraic equations §2§ and ^ are solved using the quasi- 
Newton method. Thereby we rewrite ^ and ^ in the form 

9(a) = 

where a is the vector of unknowns. Due to the fact that no more than two variables are 
coupled to each other on the left-hand side of equations ^ and ([3]), we first solve system 
(|2| for u = [tii, u 2 , u nxy ] and <p — [0i 3 4>2, •••> 0n x J i-e. the first and the third equations 
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in ^ and define 



a = (u u <f>l,U 2 , 02, ...,U nxy ,(j) nxy ) E R 



In. 



The iterative Newton method is given by 



a: 



(m+l) = a (m) _ J-l( a M) y ( a (m)) 



where the superscript denotes the iteration and Jg1 : 



:!/ is the Jacobian 



J 




Owing to the structure of the Gustafsson algorithm for the SWE, the Jacobian matrix 
is either block cyclic tridiagonal or block tridiagonal. 

J~ 1 g is solved by first applying an LU decomposition to J. Then it is computed by 
backsubstitution in two stages. First z is solved from 



In the quasi-Newton method, the computationally expensive LU decomposition is per- 
formed only once every M — th time-step, where M is a fixed integer. 

Because the backsubstitution is a fast operation, the quasi-Newton method is compu- 
tationally efficient especially when the number of nonlinear iterations at each time step is 
small. Gustafsson proved in [1] that even one quasi-Newton iteration is sufficient at each 
time step. 



Lz = g, 



and then J 



g is obtained from 



u{r l g ) 



= z. 
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The quasi-Newton formula is 



a (m+i) = a M _ J- 1 ( a M)^( a M) ) w here 
J= J(a (0) ) + O(At). 

The method works when M, the number of time-steps between successive updating 
of the LU decomposition of the Jacobian matrix J, is a relatively small number. For our 
numerical experiments we took M = 6. 

The second part of the system (J2|, the second equation in ^ is solved for v = 
[vi, t> 2 , u n ] by employing the same quasi-Newton method. Thus a is defined as 

a = ( Vl ,v 2 ,...,v nxy ) G M n ^. 

In order to obtain the SWE numerical solution at t n+1 we applied the same quasi- 
Newton technique for system (J3|. This time the coupled variables were 

a = (vi, 0i, v 3 , 02, v Hxy , <j) nxy ) G M 2n - 

for the second and third equation in (J3|, while u was solved from the remaining equation. 

3 The POD version of SWE model 

Proper orthogonal decomposition provides a technique for deriving low order model of 
dynamical systems. It can be thought of as a Galerkin approximation in the spatial 
variable built from functions corresponding to the solution of the physical system at 
specified time instances. Noack et al. [17] proposed a system reduction strategy for 
Galerkin models of fluid flows leading to dynamic models of lower order based on a 
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partition in slow, dominant and fast modes. Let us denote by Y = [u 1 , u 2 ,..., u NT ] G 
W lxvXNT an ensemble of NT— time instances of the numerical solution obtained from 
ADI FD SWE scheme at t\,t2 y -,t^T for the horizontal component of the velocity. Due 
to possible linear dependence, the snapshots themselves are not appropriate as a basis. 
Instead three methods can be employed, singular value decomposition (SVD) for Y G 
W lxyXNT , eigenvalues decomposition for YY T G ^ n *y xn xy Q r eigenvalue decomposition for 
Y T Y G ^ NTxNT (see [18] and [19]) and the leading generalized eigenfunctions are chosen 
as a basis, referred to as POD basis. Error estimates for proper orthogonal decomposition 
models for nonlinear dynamical systems may be found in [2T?] . 

Here we built the POD decomposition of each variables separately and we present only 
the construction of the POD basis corresponding to u since we applied a similar procedure 
to determine the POD bases for v and 0. Taking into account that NT <C n xy , we choose 
to construct the POD basis U G M. nxyXk , k G N* by solving the eigenvalue problem 

Y T Yui = Xm, % = 1,2,.., NT, 

and retaining the set of right singular vectors of Y corresponding to the k largest singular 
values, i.e. U = fuj}f =1 , Ui = -k=Yui. 

Similarly, let V, $ G M. nxvXk be the POD basis matrices of the vertical component of 
the velocity and geopotential, respectively. Now we can approximate u, v and <fi as follows 

u(t n ) w Uu(t n ), v(t n ) w Vv(t n ), 0(t n ) w $0(t„), 

u(t n ), v(t n ), 0(t n ) G n = 1, 2, .., NT. 

The POD reduced-order system is constructed by applying the Galerkin projection 
method to ADI FD SWE discrete model ^ and ^ by first replacing u, v, <fi with their 
approximations Uu, Vv, $0, respectively, and then premultiplying the corresponding 
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equations by U T , V T and $ T . 

The resulting POD reduced system for the first step t n+ i of the ADI FD SWE scheme 

is 



u(t n+ + ^U T F 11 (u(t n+ i),4>(t n+ i)^j = u(t n ) - ^U T F 12 (u(t n ),v(tn) 

+ ^U T (\f 2 f^ T *Vv(t n ] 

N x 

V(t n+h ) + ^-V T F 21 (tv(t n+k ), V(t n+ i)) + Y VT (M^3 T * U ^n+\) 

N x 

= v(4)-^V T F 22 (v(t„),0(t n )^, 
4>{t n+i ) + ^ T F 31 U(t nH ),4,(t n+h ) \ = 4>(t n ) - ^ T F 32 U(t n ),^(t n ) ), 



(4) 



where F n , F 12 , F 21 , F 22 , F 31 , F 32 ifxl^r are defined by 



F n (u, 0) = (Uu) * (AJJu) + -($0) * (Agt), 
F 12 (u,v) = (Vv) * (^C/u),F 21 (u,v) = (C/u) * (A x Vv), 



F 22 (v» = (Fv) * (A y Vv) + -($0) * (A/l>0), 



(5) 



F 31 (u, cj>) = -($0) * (AxUu) + (C/u) * {A£4>), 
F 32 (v,4>) = l -{®4>) * (A y Vv) + (Vv) * (A y $4>). 



The second step of the POD reduced system for the ADI FD SWE scheme is depicted 
below 
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^u T hi(v(t n+l )^(t n+ ^ 

v(Wi) + ^V r i^v(Wi),0(Wi)J = v(t n+| ) - ^V r i^u(t„ + i),v(t n+ i)J (6) 

-^V T ([f 1 f_f] T *f/a(t n+| )^ 

N x 

fctn+l) + Y $T ^ 2 (v(t n+1 ),0(W 1 )^) = 0(i n+ i) - ^$ T F 31 ^(t n+ i),0(t n+ i)^). 

The initial conditions are obtain by multiplying the following three equations with 
U T , V T , $ T 

u(f x ) « tfu(*i), v(ti) « Fv^), « $0(ti). 



We get 



ufa) « C/ T u(t!), v(*i) « V T v(*i), « $ T 0(tj 



Next we define Ai,A 2 G IR"^** 1 such as 



AiO,*) = * V(:,i), A 2 (:,i) = [fj,^) 1 * U(:,i) 



■ — 1, .., k 



-v 

iv. 



and the linear terms in (4) and (6), U T (\ f, /^.., / ] T * Vv^ and y T f [ /, /^.., / ] T * Uuj 

N x N x 



can be rewritten as U T A\ v and A 2 u respectively. 

The coefficient matrices U T A\, V T A 2 e M fcxfc defined in the linear terms of the POD 
reduced system as well as the coefficient matrices in the nonlinear functions (i.e. A X U, A y U, 
A X V, A y V, A x § , A y § e IR nxfc grouped by the curly braces in ([5])) can be precomputed, 
saved and re-used in all time steps of the interval of integration [0, i/]. However, per- 
forming the componentwise multiplications in ^ and computing the projected nonlinear 
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terms in Q and (J6]) 




U?F u (u,4>), f/ T F 12 (G,v), V T F 21 (u,v) 



(7) 



V T F 22 (v,4>), $ T F 31 (Q,0), $ T F 32 (v,0) 



still have computational complexities depending on the dimension n xy of the original 
system from both evaluating the nonlinear functions and performing matrix multiplica- 
tions to project on POD bases. If we denote the complexity for evaluating the nonlinear 



0(a(n xy ) + 4:n xy k). 

By employing the discrete empirical interpolation method we aim to remove this depen- 
dency and regain the full model reduction expected from the POD model. The projected 
nonlinear functions can be approximated by DEIM in a form that enables precomputation 
so that the computational cost is decreased and independent of the original system. Only 
a few entries of the nonlinear term corresponding to the specially selected interpolation 
indices from DEIM algorithm described in the next section must be evaluated at each 
time step. 

DEIM approximation is applied to each of the nonlinear functions F n , F 12 , F 2 ±, F 22 , F 3 i, 



F 32 defined in (JsJ) . 

4 The POD/DEIM method and its application to the 



4.1 Discrete Empirical Interpolation Method 

DEIM is a discrete variation of the Empirical Interpolation method (EIM) proposed by 
Barrault et al. [T2]. The application was suggested and analyzed by Chaturantabut and 



function Fn by a(n xy ), then the complexity for computing U T F] 



11 (u, <fi) is approximately 



ADI/SWE model 
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Sorensen in [T3], |14j . [15] . In [13], authors present an error estimate of the POD/DEIM 
method. This discrete empirical interpolation method provides an efficient way to ap- 
proximate nonlinear functions. It was also incorporated into the reduced-basis techniques 
to provide a better reduced-basis treatment (in terms of CPU time) of nonaffine and 
non-linear parameterized PDEs. DEIM was succesfully applied in conjunction with POD 
for models governing the voltage dynamics of neurons in [2T] , the integrated circuits with 
semiconductors with modified nodal analysis and drift diffusion (see [22]) and dynamics 
of the concentration of lithium ions in lithium ion batteries in [23]. In order to improve 
the stability of POD/DEIM reduced order schemes in case of a nonlinear transmission 
line, a micromachined switch and a nonlinear thermal model for a RF amplifier a few 
modifications to the DEIM based model reduction were proposed by Hochman et al. in 

Next we describe the DEIM approximation procedure applied to a nonlinear function. 
Let /:/)—>• IR n , D C IR n be a nonlinear function. If U = [ui, .., u m ), Ui G W 1 , i — 1, .., m 
is a linearly independent set, for m < n, then for r G D, the DEIM approximation of 
order m for /(r) in the space spanned by {ui} 1 ^ =l is given by 

/(r) « Uc{t), U G R nxm , c(r) G R m . (8) 

The basis U can be constructed effectively by applying the POD method on the nonlinear 
snapshots /(t**), r ti G D ( r may be a function defined from [0, T] — >■ -D, and r* 1 is the 
value of r evaluated at ti), i = 1, ..,n s , n s > 0. Next, interpolation is used to determine 
the coefficient vector c(r) by selecting m rows pi, ..,p m , Pi G N*, of the overdetermined 
linear system ^ to form a m— by— m linear system 

P T Uc(r)=P T f(r), 
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where P = [e pi , .., e p J G R nxm } e Pi = [0, -0,1,, 0, .., 0] T G R n . The DEIM approxima- 

pi 

tion of / G W 1 becomes 

f(r) « U(P T U)- l P T f(r). 

Now the only unknowns that need to be specified are the indices pi, P2, p m or the matrix 
P whose dimensions are nxm. These are determined by the following pseudo - algorithm 

DEIM: Algorithm for Interpolation Indices 

INPUT: {ui}^ C W 1 (linearly independent): 
OUTPUT: p= [px, .., p m ] G N m 

(Al). [\if)\ pi] = max\ui\,ip G M. and p\ is the component position of the largest absolute 
value of Mi, with the smallest index taken in case of a tie. 

(A2). U = [ui] G R n , P = [e Pl ] el", p= [px] G N. 

(A3). For I = 2, ..,m do 

(a) . Solve (P T U)c = P T u t for cG» w ; U, P £ W 1 *^. 

(b) . r — u\ — Uc, r G E n . 

(c) . [IV'I P/] = maa;{|r|}. 

P 



(d) . [/^[C/ u,], P^[P e pi ], p^ 



end For. 



Pi 



The DEIM procedure inductively constructs a set of indices from a linearly indepen- 
dent set. An error analysis in [TB] shows that the POD basis is a suitable choice for this 
algorithm and the order of the input basis {ui}^ C M. n according to the dominant sin- 
gular values must be utilized. Initially the algorithm searches for the largest value of the 
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first POD basis |wi| and the corresponding index represents the first DEIM interpolation 
index p\ G {1, 2, ..,n}. The remaining interpolation indices pi, I = 2, 3..,m are selected so 
that each of them corresponds to the entry of the largest magnitude of |r| defined in step 
(A3) — b. The vector r can be viewed as the residual or the error between the input basis 
ui, I = 2,3..,m and its approximation Uc from interpolating the basis {u\,U2, ..,ui-\} 
at the indices pi, p 2 , Pi-i- The linear independence of the input basis {ui} 1 £ =l guaran- 
tees that, in each iteration, r is a nonzero vector and the output indices are not 
repeating. 

An error bound for the DEIM approximation is provided in Chaturantabut and Sorensen 
[13] and [16J. An example of DEIM approximation of a highly nonlinear function defined 
on a discrete ID spatial domain can be found in [TJ], underlying the DEIM efficiency 

4.2 The DEIM SWE model 

The DEIM approximation presented earlier in this section is used to approximate the 
nonlinear terms in POD ADI SWE model described in ^ so that the nonlinear approx- 
imations will have a computational complexity proportional to the number of reduced 
variables obtained with POD. Thus, the application of DEIM in POD framework will al- 
low the construction of faster reduced order models increasing the performances of reduced 
order hierarchy models such as the one presented in Noack et al. [25] . 

Let U Fl1 G M. nxvXrn , m < n xy , be the POD basis matrix of rank m for snapshots of the 
nonlinear function Fn (obtained from ADI FD SWE scheme). Using the DEIM algorithm 
we select a set of m DEIM indices corresponding to U Fl1 , denoting by [pf 11 , .., p^Y e 
and determine the matrix Pp u G W nxyXm . The DEIM approximation of F u assumes the 
form 

so the projected nonlinear term U T Fn(u, <p) in the POD reduced system can be approx- 
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imated as 

ElgR fexm mxl 

where F™(u, $) = Pj u Fn(u, 0) and ^ G M fcxm 

Since Fn is a pointwise function (introduced in (JsJ) ) , F™ : IR fe x IR fc — > M. m can be 
defined as 

If we denote by U Fl2 , U F21 ,U F22 ,U Fsi ,U F32 G ]R n ^ xm the POD bases matrices of rank 
m for the snapshots of the nonlinear functions F 12 , F 2 i, F22, F31, F32, we obtain in a similar 
manner the DEIM approximations for the rest of the projected nonlinear terms in 



U T F l2 


(u,v) 


« u T u Fi2 {p£ i2 u Fi2 y 








s v- 

E 2 eR' sxm 


mxl 


V T F 2l 


(u,v) 


« v T u F2i (p^ 2i u F2i y 

V 


-' S v ' 






_E 3 eiR' c><m 


mxl 


V T F 22 { 




« v T u F22 (p^ 22 u F22 y 


^ v ' 






S V 


mxl 






~ <s> T u F3i (p% 3l u F3i y 


^ V v ' 






" v— 

_B 5 eiR fcxm 


mxl 


$ T F 32 ( 




V v 


_1 

-" V v ' 



£ 6 eiR' cxm mxl 
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where E 2 , E 3 , E A , E 5 , E 6 e R kxm and 

F-(u, v) = (Pj 12 Vv) * (P^AyUu), F-(u, v) = (P^Un) * (Pg^v), 
i>(v, $) = (Pj 22 Vv) * (P^AyVv) + \{PlM * (PgA* $>), 
P 3 7(u, 0) = (Pj 31 *0) * (P|AC/u) + (Pj 31 C/u) * (P|A$ 0), 

Each of the k x m coefficient matrices grouped by the curly brackets in ([9]), as well 
as Ei, i = 1, 2, .., 6 can be precomputed and reused at all time steps, so that the compu- 
tational complexity for each of the approximate nonlinear terms is 0(a(m) + 4mfc), thus 
not depending on the full-order dimension n xy . Finally, the POD/DEIM reduced system 
for the first step of ADI FD SWE model is of the form 

u(t n+| ) + ^iF-fu^:),^^^ = u(t n ) - ^E 2 F™(u(t n )Mtn)) + ^VAv(U 
v(t n+| ) + ^P 3 P 2 7(u(t n+ x),v(t n+ i)^ + ^V T A 2 u(t n+ i) = v(t n ) - Y E ^fo»)'WS) 
4>(t n+ i) + ^P 5 i^u(t n+| ),0(t n+| )) = 4>(t n ) - ^E 6 F^U(t n )^(t n )\ n = 1,..,NT-1 

(10) 

while the second step is introduced below 



fi(t nH 
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E 2 F™[ u(t n+1 ),v(t n+1 ) 



At 



v(Wi) + ~yE 4 f^[ v(t„ +1 ),0(t n+1 ; 

At 

0(t n+1 ) + — E 6 F™[ v(Ui),0(Wi) 



At. 



tf T Aiv(t n+ i) = G(t n+ i 



At 



2 



E X F£ u(t rt+ i),0(t n+ i 



0(t 



At 

T 

At 



E 3 F^[n(t n+h ),v(t n+h : 



At 



^ T A 2 u(t r 



P 5 P 3 ™ u(t n+| ), 0(t n+ O , n = 1, .., NT - 1. 



'IF 
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The initial conditions remain the same as in the case of POD reduced system. The 
nonlinear algebraic systems (10) and (10) as well as Q and ^ obtained by employing 
POD/DEIM and POD methods on ADI FD SWE were first splitted into subsystems 
according to the left-hand side of the equations where no more than two variables are 
coupled to each other and this was done in the same manner as in the case of Gustafsson 
ADI FD SWE nonlinear systems. The derived systems were solved using the Newton 
method. 



5 Numerical experiments 

In this section, we present two main experiments for the two - dimensional shallow water 
equations model to validate the fesability and efficiency of the POD/DEIM method in 
comparison with POD technique. For all tests we derived the initial conditions from the 
initial height condition No. 1 of Grammeltvedt [26] i.e. 

h(x, y) = H + H 1 + tanh( 9 ^^) + tf 2 sech 2 (s^^) ^(^) , 

< x < L, 0<y <D. 

The initial velocity fields were derived from the initial height field using the geostrophic 
relationship 

f —9\dh ( g\dh 

u= {t)w v= {7)s- x - 

Figure [T] depicts the initial geopotential isolines and the geostrophic wind field. 

The constants used were L = 6000/cm, D = 4400fcm, / = 10 _4 s _1 ,/3 = 1.5 ■ 
10~ 11 s~ 1 m~ 1 , g = 10ms~ 2 , H = 2000m, Hi = 220m, H 2 = 133m. 

For the first test, the domain was discretized using a mesh of 301 x 221 points, with 
Ax = Ay = 20km. Thus the dimension of the full-order discretized model is 66521. The 
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Figure 1: Initial condition: Geopotential height field for the Grammeltvedt initial condi- 
tion (left). Wind field (the velocity unit is lkm/s) calculated from the geopotential field 
by using the geostrophic approximation (right). 

integration time window was 24/t and we used 91 time steps (NT = 91) with At = 960s. 

ADI FD SWE scheme proposed by Gustafsson in p] was first employed in order to 
obtain the numerical solution of the SWE model. The implicit scheme allowed us to 
integrate in time using a larger time step deduced from the following Courant-Friedrichs- 
Levy (CFL) condition 

The nonlinear algebraic systems of ADI FD SWE scheme were solved with the Quasi- 
Newton method and the LU decomposition was performed only once every 6 — th time 
step. The SWE solutions at t — 24h are illustrated in Figure [2] 

The POD basis functions were constructed using 91 snapshots obtained from the 
numerical solution of the full - order ADI FD SWE model at equally spaced time steps 
in the interval [0, 24/t]. Figure [3] shows the decay around the eigenvalues of the snapshot 
solutions for u, v, <p and the nonlinear snapshots Fn, F\ 2) F 2 i, F 22 , F 3 i, F 32 . 
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Figure 2: The geopotential field (left) and the wind field (the velocity unit is lkm/s) at 
t = t f = 24h obtained using the ADI FD SWE scheme for At = 960s. 




Figure 3: The decay around the singular values of the snapshots solutions for u, v, (J) and 
nonlinear functions for At = 960s. 

The dimension of the POD bases for each variable was taken to be 35, capturing 
more than 99.9% of the system energy. We applied the DEIM algorithm for interpolation 
indices to improve the efficiency of the POD approximation and to achieve a complexity 
reduction of the nonlinear terms with a complexity proportional to the number of reduced 
variables. Figure [4] illustrates the distribution of the first 40 spatial points selected from 
the DEIM algorithm using the POD bases of F 3 i and F 32 as inputs. 
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Figure 4: First 40 points selected by DEIM for the nonlinear functions F$i (left) and F32 
(right) 

We emphasize the performances of POD/DEIM method in comparison with the POD 
approach using the numerical solution of the ADI FD SWE model. Figure [5] depicts the 
grid point local error behaviors between POD, POD/DEIM ADI SWE solutions and ADI 
FD SWE solutions, where we used 90 DEIM points. 

Using the following norms 



1 X \\w ADI FD (:,t)-w POD ADI (:,t)\\ 2 



NT 



E 

i=l 



1 ^ \ \ W ADI FD(.jj _ w POD/DEIM ADI^^ 



\ W ADI ™( : ,i)|| 2 



' NT 



E 

i=l 



\ W ADI FDr : i )U 



i = 1,2, ..,tf we calculated the average relative errors in Euclidian norm for all three 
variables of SWE model w = u,v,4>. The results are presented in Table [TJ 





POD ADI SWE 


POD/DEIM ADI SWE 




7.127e-005 


1.106e-004 


E u 


4.905e-003 


6.189e-003 


E„ 


6.356e-003 


9.183e-003 



Table 1: Average relative errors for each of the model variables. The POD bases dimen- 
sions were taken 35 capturing more than 99.9% of the system energy. 90 DEIM points 
were chosen. 
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Figure 5: Local errors between POD, POD/DEIM ADI SWE solutions and the ADI FD 
SWE solutions at t = 24h (At = 960s). The number of DEIM points was taken 90. 

In addition to the ADI FD SWE scheme we propose an Euler explicit FD SWE scheme 
as the starting point for a POD, POD/DEIM reduced model. The POD bases were 
constructed using the same 91 snapshots as in the POD ADI SWE case, only this time 
the Galerkin projection was applied to the Euler FD SWE model. The DEIM algorithm 
was used again and the numerical results are provided in Table |2j This time we employed 
the root mean square error calculation in order to compare the POD and POD/DEIM 
techniques at time t = 24h. 





ADI SWE 


POD ADI SWE 


POD/DEIM ADI SWE 


POD EE SWE 


POD/DEIM EE SWE 


CPU time seconds 


73.081 


43.021 


0.582 


43.921 


0.639 


RMSEj, 




5.416e-005 


9.668e-005 


1.545e-004 


1.792e-004 


RMSE U 




1.650e-004 


2.579e-004 


1.918e-004 


3.126e-004 


RMSE V 




8.795e-005 


1.604e-004 


1.667e-004 


2.237e-004 



Table 2: CPU time gains and the root mean square errors for each of the model variables 
at t — tf. The POD bases dimensions were taken as 35 capturing more than 99.9% of the 
system energy. 90 DEIM points were chosen. 

Applying DEIM method to POD ADI SWE model we reduced the computational time 
by a factor of 73.91. In the case of the explicit scheme the DEIM algorithm decreased the 
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CPU time by a factor of 68.733. The POD/DEIM EE SWE model was solved using the 
Runge-Kutta-Fehlberg method (RKF45). Due to the large number of spatial discretiza- 
tion points, small number of time steps and only one Newton iteration threshold imposed 
when solving the nonlinear algebraic systems of POD and POD/DEIM ADI SWE schemes 
made these two implicit schemes faster than than the POD and POD/DEIM EE SWE 
explicit schemes. The numerical results obtained showed also that the implicit schemes 
are slightly more accurate than the explicit ones. 

Figure [6] illustrates the efficiency of POD/DEIM methods as a function of spatial 
discretization points. Once the number of spatial discrete points is larger than 10000 the 
POD/DEIM schemes are faster than the POD schemes by a factor of 10, for 90 points 
selected by DEIM algorithm. 

CPU time vs. number of spatial discretization points 
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Figure 6: Cpu time vs. Spatial discretization points; POD DIM = 35, No. DEIM points 
= 90. 

Next we carried out a second experiment to test the performances of POD/DEIM 
methods. We increased the number of time steps as well as the number of snapshots 
used to generate the POD bases. Thus, we took NT = 181 and the number of snapshots 
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ADI SWE 

POD ADI SWE 

POD/DEIM ADI SWE 



POD EE SWE 
POD/DEIM EE SWE 




n s = 181. Due to the rather demanding memory requirement we had to decrease the 
number of spatial discretization points. As a consequence we choose N x = 151 and 
N y — 111, with Ax = Ay = 40km and At = 480s. 

We solved again the SWE model using the ADI FD SWE scheme in order to generate 
the 181 snapshots required for POD and POD/DEIM reduced systems. This time the 
Courant-Friedrichs-Levy (CFL) condition was 

y/ghA < 1.797. 
Ax 

The results obtained are similar with the ones obtained for a CFL condition y/gh(^) < 
7.188 underlying the performance of fully implicit Gustafsson scheme. The geopotential 
and wind field at final time tf = 24h are depicted in Figure [7j 

Contour of geopotential at time t f = 24h Wind field at time t, = 24h 




x(km) x(km) 



Figure 7: The geopotential field (left) and the wind field (the velocity unit is lkm/s) at 
t = t f = 24h obtained using the ADI FD SWE scheme for At = 480s. 

Figure [8] shows the decay of the singular values of the snapshot solutions for u, v, (j) 
and the nonlinear snapshots F 12 , F 21 , F 22 , F 31 , F 32 . We noticed that the singular 
values of F$i and F32 are decreasing slowly when compared with the singular values of 
the other nonlinear functions. 

The dimension of the POD bases for each variable was taken 35. Next we apply the 
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Figure 8: Singular values of the snapshots solutions for u, v, <p an d nonlinear functions 
for At = 480s. 



DEIM algorithm using as input the POD bases corresponding to the nonlinear functions. 
The first 40 points selected by the discrete empirical interpolation method for F31 and 
F 32 are illustrated in Figure [9j 
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Figure 9: First 40 points selected by DEIM for the nonlinear functions FF31 (left) and 
FF32 (right) 



Next we determineed solutions of the POD ADI SWE model and the POD/DEIM 
ADI SWE model using 80 DEIM points. The solutions of POD/DEIM implicit scheme 



are very accurate, local errors depicted in Figure 10 average relative errors in Table [3 
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and RMSE results in Table |4| confirm it and showing that POD/DEIM ADI SWE scheme 
is much faster and almost as accurate as POD ADI SWE scheme. 

Compared with the first experiment we reduced the number of spatial discretization 
points by a factor of 4. This does not affect the magnitude of the local errors even if 
they were decreased for both POD and POD/DEIM ADI methods with factors between 3 
and 4, when compared with the results obtained in the first case. The increased number 
of time steps and snapshots is responsible for improving the accuracy of the solutions. 



This can be observed in Figure 11 where we illustrate the local errors for the POD and 
POD/DEIM ADI SWE solutions using the same configuration as in the second experiment 
but having decreased the number of time steps and snapshots at 91. 



POD errors o„ 



DEIM/POD errors 




Figure 10: Local errors between POD, POD/DEIM ADI SWE solutions and the ADI FD 
SWE solutions at t = 24h (At = 480s). The number of DEIM points was taken 80. 



Table |3j compares the accuracy of POD and POD/DEIM ADI SWE schemes measuring 
the average relative errors of the solutions with respect to the ADI FD SWE solutions. 

Once again we calculate the solution of SWE model using the POD and POD/DEIM 
EE SWE schemes. By employing the DEIM method on the POD ADI SWE model we 
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Figure 11: Local errors between POD, POD/DEIM ADI SWE solutions and the ADI FD 
SWE solutions at t — 24h (At = 960s). The number of snapshots is 91 and the number 
of DEIM points was taken 80. 





POD ADI SWE 


DEIM/POD ADI SWE 




2.648e-005 


3.073e-005 




1.279e-003 


1.292e-003 


E v 


2.207e-003 


2.471e-003 



Table 3: Average relative errors for each of the model variables at t = tf,At = 480s. The 
POD bases dimensions were taken 35 capturing more than 99.9% of the system energy. 
80 DEIM points were chosen. 

reduced the CPU time with a factor of 13. In the case of explicit scheme the DEIM 
algorithm decreased the computational time by a factor of 21. The adaptative Runge- 
Kutta-Fehlberg method involved in the explicit reduced order models (ROMs) was faster 
than the Newton method used to solve the nonlinear algebraic systems in implicit ROMs 
mostly because we doubled the number of time steps and thus the RKF45 didn't need 
to generate a large amount of intermediary time steps as it did in the first experiment 
in order to generate an accurate solution. From Table [4] we notice that the RMSEs for 
both implicit and explicit POD/DEIM schemes are almost similar with the ones generated 
by the implicit and explicit POD systems with respect to the ADI FD SWE numerical 
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solutions. 

Collecting the results obtained from experiments 1 and 2 we conclude that the POD 
and POD/DEIM ADI SWE schemes are more accurate than the POD and POD/DEIM 
EE SWE schemes. 





ADI SWE 


POD ADI SWE 


DEIM/POD ADI SWE 


POD EE SWE 


DEIM/POD EE SWE 


CPU time 


24.353 


8.848 


0.686 


8.019 


0.386 


RMSEq 




1.607e-005 


1.743e-005 


3.999e-005 


4.074e^005 


RMSE U 




4.076e-005 


4.233e-005 


5.089e-005 


5.780e-005 


RMSE V 




2.397e-005 


2.755e-005 


4.613e-005 


4.919e-005 



Table 4: CPU time gains and the root mean square errors for each of the model variables 
at t = tf, At = 480s. The POD bases dimensions were taken 35 capturing more than 
99.9% of the system energy. 80 DEIM points were chosen. 



Next we evaluate the efficiency of POD/DEIM method as a function of POD dimen- 



sion. Figure 12 gives the root mean square errors of and the corresponding average 
CPU times for different dimensions of POD and DEIM approximations. 



CPU time 



■ FULL 

- POD ADI SWE 
POD EE SWE 

■ DEIM/POD ADI SWE 50 
DEIM/POD ADI SWE 120 

- DEIM/POD EE SWE 50 
DEIM/POD EE SWE 120 




Root mean square error of 4 



- POD ADI SWE 
POD EE SWE 

■ DEIM/POD ADI SWE 50 
DEIM/POD ADI SWE 120 

■ DEIM/POD EE SWE 50 
DEIM/POD EE SWE 120 




POD dimension 



15 20 25 30 

POD dimension 



Figure 12: CPU time of the full system, POD reduced systems and POD - DEIM reduced 
systems (left); Root mean square error of calculated for POD/DEIM and POD reduced 
systems with respect to ADI FD SWE solutions 



The results in Tables |5| and |6| show the performances of POD/DEIM method. Once 
the POD dimension exceeds 15(5) in case of DEIM/POD ADI (EE) SWE scheme the 
CPU time is decreased at least by a factor of 10. When DEIM dimension reached 50, we 
notice that RMSE results between POD and POD/DEIM are almost identical. All the 
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methods performed well when POD dimension exceeded 25 leading to RMSE results of 
order 0(1(T 5 ). 



PODDIM 


P( ID ADI SWE 


DEIM/POD ADI SWE50 


DEM/PI ID ADI SWE120 


POD EE SWE 


DEIM/POD EE SWE50 


DEIM/POD EE SWE120 


1 


0.0009 


0.0003 


0.0002 


0.0004 


0.0001 


0.0002 




0.0204 


0.0109 


0.0111 


0.0235 


0.0014 


0.0014 


10 


0.1896 


0.0244 


0.0249 


0.1838 


0.0066 


0.0071 


15 


0.6478 


0.0599 


0.0615 


0.6211 


0.0238 


0.0251 


20 


1.5336 


0.1245 


0.1283 


1.4795 


0.0586 


0.0615 


25 


2.9991 


0.2326 


0.2403 


2.8901 


0.1165 


0.1237 


30 


5.1956 


0.3840 


0.4013 


5.0250 


0.2202 


0.2292 


35 


8.8480 


0.6718 


0.6971 


8.0190 


0.3773 


0.3992 



Table 5: Comparison between CPU times of POD and DEIM/POD implicit and explicit 
schemes. The computational time of the full model (ADI FD SWE) was 24.3530. 



PODDIM 


POD ADI SWE 


DEIM/POD ADI SWE50 


DEIM/POD ADI SWE120 


POD EE SWE 


DEIM/POD EE SWE50 


DEIM/POD EE SWE120 


1 


1.182c-003 


1.182c-003 


1.182c-003 


1.1826-003 


1.182c-003 


1.1826-003 




5.018c-004 


5.087c-004 


5.047e-004 


5.037o-004 


5.110c-004 


5.070o-004 


10 


4.526c-004 


5.122c-004 


4.609c-004 


4.553c-004 


5.155c-004 


4.635c-004 


15 


1.587e-004 


1.6826-004 


1.589o004 


1.584o-004 


1.669c-004 


1.595e-004 


20 


1.009c-004 


l.lOlc-004 


1.012c-004 


9.876c-005 


1.097c-004 


9.8766-005 


25 


5.075C-005 


5.526c-005 


5.438c-005 


5.9096-005 


6.319e-005 


6.281e-005 


30 


2.765c-005 


2.898o-005 


2.8006-005 


4.788o-005 


5.149c-005 


4.923o-005 


35 


1.607c-005 


1.846c-005 


1.835c-005 


3.999c-005 


4.439e-005 


4.175e-005 



Table 6: Comparison between RMSE of POD and DEIM/POD implicit and explicit 
schemes. 



Numerical experiments carried out for a 1-day integration showed that the POD/DEIM 
ADI SWE scheme as well as the POD/DEIM EE SWE discrete model conserve the average 
height of the free surface and the potential enstrophy while another integral invariant 
of the SWE model, the total energy, is not preserved largely due to the absence of a 
staggered C-grid in the numerical discretization. Arakawa in [27] showed that when the 
finite difference Jacobian expression for the advection term is restricted to a form which 
properly represents the interaction between grid points the computational instability is 
prevented thereby preserving all the integral invariants. Thus the POD and POD/DEIM 
systems behave in a similar way as the ADI FD SWE full scheme in the matter of integral 



invariants conservations and Figure 13 depicts their evolution in time. 

Tables [7j - [9] present integral invariants measures for all the schemes involved in this 
study using max — min evaluation and Euclidian norm with respect to SWE invariants 
calculated with ADI FD SWE full scheme. 
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- AD I SWE 

- POD AD I SWE 

- DEIM/POD ADI SWE " 
POD EE SWE 
DEIM'POD EE SWE 




7.8638 
7.8536 
7.8534 
7.8532 
7.853 



Potential Enstrophy 



Average Height of the Free Surface 



- ADI SWE 

- POD ADI SWE 

- DEIM/POD ADI SWE 
POD EE SWE 
DEIM/POD EE SWE 




- ADI SWE 

- POD ADI SWE 

- DEIM/POD ADI SWE 
POD EE SWE 
DEIM/POD EE SWE 




Figure 13: Shallow Water Equations invariants, POD dimension =35, No of DEIM points 
= 120 





ADI SWE 


POD ADI SWE 


DEIM/POD ADI SWE 


POD EE SWE 


DEIM/POD EE SWE 


Max-Min 


0.0017 


0.0063 


0.0125 


0.0060 


0.0237 


Norm 


0.0000 


0.0331 


0.0693 


0.0348 


0.0795 



Table 7: Average Height of the conservation of the mass 



In Figure 14 , the Pearson correlation coefficient defined below is used as an additional 



metric to evaluate the quality of POD/DEIM schemes 



where 



cov 



12 



i = \,..,NT, 



a 



J—n X y 



J—n X y 



{Wij - Wj) 2 , o-a = (w°f eme - W scheme j ) , i = l,..., NT, 



3=1 



3 — T^xy 



cov 12 = ( W i,j ~ W i) (W t s f eme - W scheme ^j , 



i = l,...,NT, 



where W = u,v, <fr represents the ADI FD SWE solution and W 



scheme „. scheme „, scheme Ascheme 



U 



V 



the solution calculated with one of the following schemes: POD ADI SWE, POD/DEIM 
ADI SWE, POD EE SWE and POD/DEIM EE SWE using 50 and 120 DEIM points. Wj 



and W scheme j are corresponding means over the simulation period [0,tf] at spatial node 
J- 
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ADI SWE 


POD ADI SWE 


DEIM/POD ADI SWE 


POD EE SWE 


DEIM/POD EE SWE 


Max-Min 


0.0012 


0.0011 


0.0014 


0.0011 


0.0008 


Norm 


0.0000 


0.0014 


0.0023 


0.0024 


0.0040 



Table 8: Potential Enstrophy 





ADI SWE 


POD ADI SWE 


DEIM/POD ADI SWE 


POD EE SWE 


DEIM/POD EE SWE 


Max-Min 


5.6440e+016 


5.7392e+016 


5.6140e+016 


6.0441e+016 


6.3439e+016 


Norm 


0.0000e+000 


1.4289e+016 


3.3597e+016 


2.9898e+016 


5.3457e+016 



Table 9: Total Energy 




Figure 14: Correlation coefficients for the SWE variables, POD dimension =35. 

We also tested the SWE discrete schemes in context of parametric variation of the 
Coriolis parameter. For these tests we used a mesh of 301 x 221 points and we run the 
ADI FD scheme to generate 91 snapshots using / = 10~ 4 . Next we derived the POD 
bases and than we solved the POD ADI SWE and POD/DEIM SWE discrete problems 



for / = 3 ■ 10 4 . Figure 15 depicts the grid point local error behaviours between POD, 
POD/DEIM ADI SWE solutions and ADI FD SWE solutions obtained by assigning the 
Coriolis parameter the value 3 • 10~ 4 . The DEIM scheme used 90 interpolations points. 



Tables 10 and 11 contain the average relative errors and RMSEs for u, v and 0. Cor- 
relation coefficients are shown in Figure [TBJ 

We can improve the accuracy of POD, POD/DEIM solutions and reduce the com- 
putational time by taking into account different POD bases dimensions according to the 
eigenvalues decay. The final numerical test presented here is dedicated to the effect of 
different dimensions of POD bases in POD and POD/DEIM Euler explicit SWE schemes. 
As we expected, after a short analysis of eigenvalues spectrum we chose the dimensions of 
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Figure 15: Local errors between POD, POD/DEIM ADI SWE solutions and the ADI FD 
SWE solutions at t — 24h (At = 960s) in context of parameter settings. The number of 
DEIM points was taken 90. 



POD bases to be 31, 32 and 34. The results are presented in Table 12 and we notice that 



the errors are smaller then the ones in Table [2] where we used 35 modes for each POD 
basis. 



6 Conclusions 

To obtain the approximate solution in case of both POD and POD/DEIM reduced sys- 
tems, one must store POD or POD/DEIM solutions of order O(kNT), k being the POD 
bases dimension and NT the number of time steps in the integration window. The coeffi- 
cient matrices that must be retained while solving the POD reduced system are of order 
of 0(k 2 ) for projected linear terms and 0(n xy k) for the nonlinear term, where n xy is the 
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POD ADI SWE 


POD/DEIM ADI SWE 




0.124508 


0.14588 


E u 


0.22403 


0.27408 


E v 


1 .4169e-003 


1.6294e-003 



Table 10: Average relative errors for each of the model variables. The POD bases dimen- 
sions were taken to be 35 capturing more than 99.9% of the system energy. 90 DEIM 
points were chosen. 





POD ADI SWE 


POD/DEIM ADI SWE 


RMSE^ 


5.6992e-004 


8.3716e-004 


RMSE U 


8.8999e-004 


1.2342e-003 


RM SE V 


3.56849e-004 


5.2638e-003 



Table 11: Root mean square errors for each of the model variables at time t = 2 Ah. The 
POD bases dimensions were taken 35 and 90 DEIM points were chosen. 

space dimension. 

In the case of solving POD/DEIM reduced system the coefficient matrices that need 
to be stored are of order of 0(k 2 ) for projected linear terms and 0(mk) for the nonlinear 
terms, where m is the number of DEIM points determined by the DEIM indexes algo- 
rithm, m <C n xy . Therefore DEIM improves the efficiency of the POD approximation and 
achieves a complexity reduction of the nonlinear term with a complexity proportional to 
the number of reduced variables. 

We proved the efficiency of DEIM using two different schemes, the ADI FD SWE fully 
implicit model and the Euler explicit FD SWE scheme. We noticed, as we expected, that 
POD/DEIM CPU time is most sensitive to the number of spatial discretization points. 
The largest reduction of the CPU time was obtained in first experiment where it was 
reduced by a factor of 73.91 when using the POD/DEIM ADI SWE scheme while in the 
case of POD/DEIM EE SWE model we decreased the CPU time by a factor of 68.733. 
Also we noticed that the approximation errors of POD/DEIM and POD reduced systems 
are almost identical once the dimension of DEIM attained the value of 50, for any of 
the methods used, either explicit or implicit. In the second experiment, we increased the 
number of time steps and snapshots and consequently the solutions accuracy was higher 
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Figure 16: Correlation coefficients for the SWE variables in parameter settings, POD 
dimension =35. 





POD EE SWE (35 modes) 


POD/DEIM EE SWE (3. r , modes) 


POD EE SWE (31.32.34 modes) 


POD/DEEM EE SWE (31,32,34 modes) 


RMSEj, 


1.545e-004 


1.792e-004 


1.1605e-004 


1.4246e-004 


RMSE U 


1.918e-004 


3.126e-004 


9.3842e-005 


L712e-004 


RAISE,, 


1.667e-004 


2.2374e-004 


9.551e-005 


1.07916-004 



Table 12: The root mean square errors for each of the model variables at t = 24h. 
Comparison between the errors calculated with the same number of modes (first and 
second columns) and different number of modes (third and fourth columns). We used a 
mesh of 301 x 221 points and 90 DEIM points were chosen. 



in comparison with the results obtained in the first experiment. 

In future research we plan to apply the DEIM technique to different inverse problems 
such as POD 4-D VAR of the limited area finite element shallow water equations and 
adaptive POD 4-D VAR applied to a finite volume SWE model on the sphere. 

We are also interested to compare the Discrete Empirical Interpolation Method with 
the 'Best Points' Interpolation (BPIM) one in proper orthogonal decomposition frame- 
work applied to SWE equation. BPIM was proposed by Nguyen et al [2S] where the 
interpolation points are defined as a solution of a least-squares minimization problem. 
Thus, BPIM replaces the greedy algorithm used in Empirical Interpolation Method(EIM) 
by an optimization problem which provides higher accuracy at the cost of greater compu- 
tational complexity. For instance Galbally et al. [29] made a comparison between gappy 
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POD, EIM and BPIM techniques for a nonlinear combustion problem governed by an 
advection diffusion PDE. 
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Appendix 

This appendix contains a definition for the componentwise multiplication operation and a 
more symbolic representation of the Gustafsson's nonlinear ADI finite difference shallow 
water equations scheme defined in (2-3), Section 2. The componentwise multiplication is 
also known as the Hadamard product. 
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Definition The componentwise multiplication is a binary operation that takes two 
matrices of the same dimensions, A, B e M. mxn , m,n e N* ; and produces another matrix 
of the same dimensions A* B G R mxn with elements satisfying the following relationship 

{A*B) iJ = (A) i j.(B) i j. 

Now, recall that the number of mesh points is n xy = N x ■ N y . We shall denote by 
Wij(i = 1,..,N X ; j = l,..,N y ] n = 1,..,NT) an approximation to w(iAx, jAy, nAt) = 
(u,v,(f))(iAx,jAy,nAt). The basic difference operators are 



D 0x w? j = (w2 ¥1J -wZ. liS )/(2Ax), 



= «- <_!,,-) Ax, 

respectively, with similar definitions for _D j/? We also define the operators Pfy 

and g^- by 

^. = At/2(A(< J .)^o. + C'j 1) ), 

(2)>| 



g?,. = At/2(S«)^ + Cf ). 



with A, _B defined at the beginning of Section (2), 



A = < 



D 0y , j = 2,..,N y -l, 

D+y, j = 1, 
D—y, j = Ny, 
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(owing to the boundary conditions in the y— direction) 
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and fj = f{jAy), j = 1, .., N y . 

Thus we rewrite the Gustafsson's nonlinear ADI difference scheme given in (pip 



i = l,.., N x - 1; j = 1, .., N y ; n = 1,2,.., NT- 1. 



(12) 
(13) 



The periodic boundary conditions in the x— direction allowed us to use only central dif- 
ferences to approximate the derivatives with respect to x and eliminated the need of 
calculating the SWE solutions for % = N x . 

The nonlinear algebraic equations do not apply to the v component for j = 1, j = N y , 



but we used the condition v™ 1 



0, i = l,..,N x -l, n= 1,2,.., NT. 
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